use strict;
use Getopt::Mixed 1.006, 'getOptions';
use LWP::Simple;
use LEO::set_comp;
use LEO::countregex_FA_seq;

sub uniquify(\@)
{
	my %hash;
	@hash{@{$_[0]}} = ();
	@{$_[0]} = keys %hash;	
	return @{$_[0]};
}	

sub printstep($$)
{
	print "\n *********\n Step @_[0]:\n *********\n\n    @_[1]\n\n";
}	

sub fc_grep($$)
{
	chomp(my $cnt = qx"grep -c \"$_[0]\" $_[1]");
	return $cnt;
}

our $opt_help = undef;
our $opt_missensesprot = 0;
our $opt_ignoremorbidmap = 0;
our $opt_omimsvariants = 0;
our $opt_keeptemp = 0;
our $opt_dirtemp = "tmp_sp";
our $opt_redownload = 0;

getOptions(	"omimsvariants missensesprot ignoremorbidmap help dirtemp=s keeptemp redownload " .
		   	"i>ignoremorbidmap o>omimsvariants m>missensesprot h>help d>dirtemp ".
			"k>keeptemp r>redownload");


print "\n";

die << "USAGE" if (defined $opt_help);
Usage:

    disease.pl
	
		[-r|redownload]        re-download databases from omim and swissprot
        [-o|omimsvariants]     re-parse variants from the OMIM database
        [-m|missensesprot]     re-parse swissProt variants with MIM codes
        [-i|ignoremorbidmap]   ignore MIM codes from the out of date morbid map
        [-d|dirtemp="tmp_sp"]  temporary directory for output
        [-k|keeptemp]          do not delete temporary directory
        [-h|help]              this message

USAGE

print "\tReparse OMIM codes associated with missense variants.\n" if ($opt_omimsvariants);
print "\tReparse SProt sequences with variants.\n" if ($opt_missensesprot);

# temporary directory
mkdir "$opt_dirtemp" unless -f "$opt_dirtemp";
print "\tTemporary files will be stored in $opt_dirtemp.\n" if $opt_keeptemp;


{

#_________________________________________________________________________________________
#
#	Get MIMs	
#_________________________________________________________________________________________
	if ($opt_redownload)
	{
		my @date;
		$date[0] = localtime((stat('D:\werke\database\omim\morbidmap'))[9]);
		$date[1] = localtime((stat('D:\werke\database\omim\genemap'))[9]);
		$date[2] = localtime((stat('D:\werke\database\omim\omim.txt'))[9]);
		$date[3] = localtime((stat('D:\werke\database\sprot\sprot.dat'))[9]);

		printstep(0, "Redownloading OMIM and swissprot databases:");
		# morbid map
			print "\tDownloading morbid map from ncbi...\n";
			my $rc = getstore('ftp://ncbi.nlm.nih.gov/repository/OMIM/morbidmap',
					'D:\werke\database\omim\morbidmap');
			die "Download failed.\n" unless is_success($rc);
			print "\t\tDownloaded morbid map.\n";


		#gene map
			print "\n\tDownloading gene map from ncbi...\n";
			my $rc = getstore('ftp://ncbi.nlm.nih.gov/repository/OMIM/genemap',
					'D:\werke\database\omim\genemap');
			die "Download failed.\n" unless is_success($rc);
			print "\t\tDownloaded genemap.\n";
		

		# omim database
			print "\n\tDownloading omim database from ncbi...\n";
			my $rc = getstore('ftp://ncbi.nlm.nih.gov/repository/OMIM/omim.txt.Z',
					'D:\werke\database\omim\omim.txt.Z');
			die "Download failed.\n" unless is_success($rc);
			print "\t\tDownloaded omim database.\n";
			
			print "\t\tunzipping omim.txt.Z...\n\t\t";
			die "unzip failed" unless (!system('gzip -daNf D:\werke\database\omim\omim.txt.Z')
										&& -f 'D:\werke\database\omim\omim.txt');
			print "\t\tOmim.txt unzipped.\n";
		
		# swissprot database
			print "\n\tDownloading swiss prot database from ncbi...\n";
			my $rc = getstore('ftp://ftp.ebi.ac.uk/pub/databases/sp_tr_nrdb/sprot.dat.gz',
					'D:\werke\database\sprot\sprot.dat.gz');
			die "Download failed.\n" unless is_success($rc);
			
			print "\t\tunzipping sprot.dat.Z...\n";
			die "unzip failed" unless (!system('gzip -dafN D:\werke\database\sprot\sprot.dat.gz')
										&& -f 'D:\werke\database\sprot\sprot.dat');
			print "\t\tsprot.dat unzipped.\n\n";
			
		# print dates
			print "File dates used:\n";
			$date[4] = localtime((stat('D:\werke\database\omim\morbidmap'))[9]);
			$date[5] = localtime((stat('D:\werke\database\omim\genemap'))[9]);
			$date[6] = localtime((stat('D:\werke\database\omim\omim.txt'))[9]);
			$date[7] = localtime((stat('D:\werke\database\sprot\sprot.dat'))[9]);
			print "\tOMIM morbidmap dates new/old:\n\t\t$date[0]\n\t\t$date[4]\n";
			print "\tOMIM genemap dates new/old:\n\t\t$date[1]\n\t\t$date[5]\n";
			print "\tOMIM database dates new/old:\n\t\t$date[2]\n\t\t$date[6]\n";
			print "\tSwissProt database dates new/old:\n\t\t$date[3]\n\t\t$date[7]\n\n";
	}
	else
	{
		my @date;
		print "File dates used:\n";
		$date[0] = localtime((stat('D:\werke\database\omim\morbidmap'))[9]);
		$date[1] = localtime((stat('D:\werke\database\omim\genemap'))[9]);
		$date[2] = localtime((stat('D:\werke\database\omim\omim.txt'))[9]);
		$date[3] = localtime((stat('D:\werke\database\sprot\sprot.dat'))[9]);
		print "\tOMIM morbidmap dates:\n\t\t$date[0]\n";
		print "\tOMIM genemap dates:\n\t\t$date[1]\n";
		print "\tOMIM database dates:\n\t\t$date[2]\n";
		print "\tSwissProt database dates:\n\t\t$date[3]\n";
	}
	


	printstep(1, "Getting list of MIM codes from the OMIM morbid map:");
	chomp (my @OMIM_morbid_MIMs = qx'get_morbid_MIM.pl ..\database\OMIM\morbidmap');


	printstep(2, "Getting list of MIM codes from the OMIM gene map:");
	chomp (my @OMIM_gene_map_morbid_MIMs = qx'get_gene_map_morbid_MIM.pl ..\database\OMIM\genemap');

	my $records = fc_grep '^*RECORD*', '..\database\OMIM\OMIM.txt';
	
	printstep(3, "Getting list of MIM codes\n\tassociated with OMIM missense variants in $records records:");
	if ($opt_omimsvariants or ! -f"$opt_dirtemp\\variants.list.fa")
	{
		`lookup_variants_in_OMIM.exe < ..\\database\\OMIM\\OMIM.txt > $opt_dirtemp\\variants.list.fa`;
	}
	else
	{
		print "\tUse OMIM variant records retrieved earlier in variants.list.fa.\n";
	}
	print "\n    Extracting OMIM codes for missense variants:\n\n";
	chomp (my @OMIM_variant_MIMs = `Perl -S extract_regex.pl -s=">(\\d+)" $opt_dirtemp\\variants.list.fa`);



#_________________________________________________________________________________________
#
#	Merge and save MIMs
#_________________________________________________________________________________________
	printstep(4, "Merging the sets of MIM codes:");
	my ($common, $morbid_unique, $genemap_unique, $variant_unique) = 
											SetCompare3(	@OMIM_morbid_MIMs, 
															@OMIM_gene_map_morbid_MIMs, 
															@OMIM_variant_MIMs);
	my (@all_mims) = 	(	@OMIM_morbid_MIMs, 
							@OMIM_gene_map_morbid_MIMs, 
							@OMIM_variant_MIMs);
								

	push (@all_mims, @OMIM_morbid_MIMs) unless ($opt_ignoremorbidmap);
	@all_mims = uniquify (@all_mims);
	

	my $cmmon_cnt = @all_mims -@$morbid_unique- @$genemap_unique-@$variant_unique-@$common;
	print <<"    OUTPUT";
	@{[scalar @$morbid_unique]}	unique MIM codes in the OMIM morbid map.
	@{[scalar @$genemap_unique]}	unique MIM codes in the OMIM gene map.
	@{[scalar @$variant_unique]}	unique MIM codes with OMIM variants.
	@{[scalar @$common]}	       MIM codes common to all three.
	$cmmon_cnt	       other MIM codes in common.
      @{['_'x8]}
	@{[scalar @all_mims]}\tMIM codes altogether.
    OUTPUT

	print "\t\t(MIM codes unique to OMIM morbid map were discarded.\n" if ($opt_ignoremorbidmap);


	print("\n\tSaving MIM codes to 'all_MIMS.list':");
	open FILE_ALL_MIMS, ">$opt_dirtemp\\all_MIMS.list";
	print FILE_ALL_MIMS join ("\n", @all_mims);
	close FILE_ALL_MIMS;

#_________________________________________________________________________________________
#
#	Parse SwissProt for MIM codes and variants [optional/slow]
#_________________________________________________________________________________________
	if ($opt_missensesprot or ! -f 'database\SProt_with_MIMs.fa')
	{
		printstep(5, "Getting SwissProt sequences with associated MIM codes and variants:");
		qx'SProt_to_FA_if_MIM.pl -e2000 ..\database\sprot\sprot.dat > database\SProt_with_MIMs.fa';
	}
	else
	{
		printstep(5, "Using SwissProt sequences with associated MIM codes retreived earlier.");
	}
	printstep(6, "Extracting SwissProt sequences with disease MIM codes:");
	system('perl -S extract_seq_from_FA.pl -s="MIM\|(\d+)" database\SProt_with_MIMs.fa'.
			" < $opt_dirtemp\\all_MIMS.list > disease.fa");
	my $cnt_sequences = fc_grep '^>', "disease.fa";
	print "\t$cnt_sequences\tsequences were retrieved.\n";
	
	
#_________________________________________________________________________________________
#
#	Get variants records
#_________________________________________________________________________________________
	printstep(7, "Extracting variant records sorted by SwissProt accessions:");
	system('perl SPROT_VARIANTs_from_FA.pl '.
			" < disease.fa > $opt_dirtemp\\SProt_variants.list.fa");

	my $cnt_varrec = fc_grep '.*', "$opt_dirtemp\\SProt_variants.list.fa";
	my $cnt_varseq = fc_grep '^>', "$opt_dirtemp\\SProt_variants.list.fa";
	$cnt_varrec -= $cnt_varseq;

	print << "    OUTPUT";
	$cnt_varseq	sequences with variant records were retrieved.
	$cnt_varrec	variant records were retrieved.
    OUTPUT

#_________________________________________________________________________________________
#
#	Get which MIM codes succeeded
#_________________________________________________________________________________________
	printstep(8, "Getting the list of MIM codes used to retrieve the variant records:");
	
	my @variants = @OMIM_variant_MIMs;
#	SetCompare(@variants, \@)
#	return (\@arr_common, \@arr_unique1, \@arr_unique2);
#}	

	my ($cnt_seq1, $except_SP_1)
				= countregex_seq(@variants, 'MIM\|(\d+)', "disease.fa", 'SP\|[^|]*\|');
	my ($cnt_seq2, $except_SP_2)
				= countregex_seq_except(@OMIM_gene_map_morbid_MIMs, @$except_SP_1, 
											'MIM\|(\d+)', "disease.fa", 'SP\|[^|]*\|');
	push @$except_SP_1, @$except_SP_2;
	my ($cnt_seq3)
				= countregex_seq_except(@OMIM_morbid_MIMs, @$except_SP_1, 
											'MIM\|(\d+)', "disease.fa", 'SP\|[^|]*\|');

print <<"OUTPUT";
	$cnt_seq1	sequences retrieved using the variants MIM codes.
	$cnt_seq2	further sequences retrieved using the genemap MIM codes.
	$cnt_seq3	further sequences retrieved using the morbid map MIM codes.
    _________
	@{[$cnt_seq1 + $cnt_seq2 + $cnt_seq3]}	sequences in total.

OUTPUT



#_________________________________________________________________________________________
#
#	Check if variants are in sequences
#_________________________________________________________________________________________
	printstep(9, "Check if variant records map correctly to the sequences:");
	system ('Perl check_sprot_variants_against_sequence.pl '	.
					"disease.fa "					.
					"$opt_dirtemp\\SProt_variants.list.fa "		.
					"$opt_dirtemp\\SProt_variants.list.fa > nul");


#_________________________________________________________________________________________
#
#	parse domains from hmmpfam output
#_________________________________________________________________________________________
	printstep(10, "Calculating lengths of sequences retreived from SwissProt:");
	system('perl -S get_seq_len_FA.pl -a="SP\|([^|]+)" disease.fa '.
						"> $opt_dirtemp\\disease_seq_len.list");

	printstep(11, "HMMPfam:");
    print << "    HMMPFAM";
	disease.fa is in the base directory
	Run:
		hmmpfam -Z373408  all.hmm disease.fa > smart_pfam.hmmout
		hmmpfam -Z373408  smart.hmm disease.fa > smart.hmmout
		
	Copy:
		smart_pfam.hmmout	to database\\smart_pfam.hmmout
		smart.hmmout		to database\\smart.hmmout

	Press return when ready:
    HMMPFAM
	<STDIN>;
	

	printstep(12, "Parse domains from hmmpfam output:");
	system( 'parse_domains_from_hmmpfam '.
					'database\smart_pfam.hmmout '.			# smart + pfam output
					'database\smart.hmmout '.				# smart output
					"$opt_dirtemp\\disease_seq_len.list ".	# sequence lengths
					'0.1 '.									# e value cut of point
						"> smart_pfam_domains.list.fa");
						
						
	printstep(13, "Parse domains from hmmpfam output for SMART domains only:");
	system( 'parse_domains_from_hmmpfam '.
					'database\smart.hmmout '.				# smart output
					'database\smart.hmmout '.				# smart output
					"$opt_dirtemp\\disease_seq_len.list ".	# sequence lengths
					'0.1 '.									# e value cut of point
						"> nul");

	printstep(14, "Map variants onto domains");
	system('perl variants_in_domains.pl '.
				"$opt_dirtemp\\domains.list.fa ".
				"smart_pfam_domains.list.fa ".
				"< $opt_dirtemp\\SProt_variants.list.fa ".
					"> $opt_dirtemp\\SProt_variants_domains.list.fa");

	printstep(15, "Calculate domain statistics:");

	my $missed_variants = fc_grep 'empty_', "$opt_dirtemp\\SProt_variants_domains.list.fa";
	print "\t$cnt_varrec\tmissense mutations in total\n";
	print "\t$missed_variants\tmissense mutations did not map onto any domain\n";
	printf "\t%.2f%% of missense mutations mapped onto a predicted domain.\n",
			 ($cnt_varrec - $missed_variants) * 100 / $cnt_varrec;

	printstep(16, "Mapping diseases to domains for smart:");
	system ('perl sp_omim_by_domain.pl '.
				"$opt_dirtemp\\all_mims.list ".
				"disease.fa ".
				'..\database\omim\genemap '.
				'..\database\omim\morbidmap '.
				'..\database\omim\omim.txt '.
				"$opt_dirtemp\\SProt_variants_domains.list.fa".
				">SMART_disease_domains.list");

	print "\n\n";
	print "_"x79;
	print "\n\n\tOutput Files:\n";
	print "_"x79;
	print "\n\n\n";
	print "\tSMART_PFam domain/mutations \tin \"smart_pfam_domains.list.fa\".\n";
	print "\tSMART domain/diseases \t\tin \"SMART_disease_domains.list\".\n";
	print "\tSequences linked with diseases \tin \"disease.fa\".\n";

	# remove temporary directory
	unless ($opt_keeptemp)
	{
		unlink glob("$opt_dirtemp\\*") or die;
		rmdir "$opt_dirtemp";
	}
	

}	

																									

